//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "ephi.hpp"

// look at sphere of charge e-field with non-uniform charge distribution rho=k*(1-(r/R)**2)**2

class myElement : public StaticElement
{
public:
  myElement (const vect3d& pos, prec_t R, prec_t Q) : pos(pos), R(R), Q(Q) { }

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield)
  {
    vect3d r = pos - this->pos;
    prec_t r2 = r*r;
    prec_t rlen = prec_t_sqrt(r2);
    prec_t r3 = r2 * rlen;

    if ( rlen < 1e-10 )
      {
        return;
      }
    else if ( rlen >= R )
      {
        efield += (Q / r3) * r;
      }
    else
      {
        prec_t r5 = r3 * r2;
        prec_t r7 = r5 * r2;
        prec_t R2 = R * R;
        prec_t R4 = R2 * R2;
        prec_t R5 = R4 * R;
        prec_t R7 = R5 * R2;
        prec_t R9 = R7 * R2;
        efield += (Q * (35 * r3 / (8 * R5) - 21 * r5 / (4 * R7) + 15 * r7 / (8*R9))) * r.normal();
      }
  }

  void addPotential (const vect3d& pos, prec_t& potential)
  {
    vect3d r = pos - this->pos;
    prec_t r2 = r*r;
    prec_t rlen = prec_t_sqrt(r2);
    prec_t r4 = r2 * r2;
    prec_t r6 = r4 * r2;
    prec_t r8 = r6 * r2;
    if ( rlen >= R )
      {
        potential += Q / rlen;
      }
    else
      {
        prec_t R2 = R * R;
        prec_t R4 = R2 * R2;
        prec_t R5 = R4 * R;
        prec_t R7 = R5 * R2;
        prec_t R9 = R7 * R2;

        potential += -Q * (35 * r4 / (32*R5) - 21 * r6 / (24 * R7) + 15 * r8 / (64 * R9) - 93 / (64 * R));
      }
  }

  void draw(Screen& screen, const Screen::color& c)
  {
  }

  std::string get_povray_script ()
  {
    return std::string();
  }

  bool inContact (const vect3d& pos)
  {
    return false;
  }

private:
  const vect3d& pos;
  prec_t R;
  prec_t Q;
};

class MyBuilder : public OctreeBuilder
{
public:
  MyBuilder(const vect3d& p0, prec_t R, prec_t q) :
      OctreeBuilder (vect3d(-1, -1, -1), vect3d(1, 1, 1), 1, 1, 100, 0.25),
      p0(p0), R(R), q(q)
  {
    prec_t R2 = R * R;
    prec_t R5 = R2 * R2 * R;
    prec_t R7 = R5 * R2;
    prec_t R9 = R7 * R2;

    f35_8R5 = 35 / (8 * R5);
    f21_4R7 = 21 / (4 * R7);
    f15_8R9 = 15 / (8 * R9);
    f35_4R5 = 35 / (4 * R5);
    f21_R7 = 21 / R7;
    f45_4R9 = 45 / (4 * R9);
    f42_R7 = 42 / R7;
    f45_R9 = 45 / R9;
    f90_R9 = 90 / R9;
    f35_32R5 = 35 / (32 * R5);
    f21_24R7 = 21 / (24 * R7);
    f15_64R9 = 15 / (64 * R9);
    f93_64R = 93 / (64 * R);
  }

  void addValues (const vect3d& pos, CubicScalar *scalars, CubicVector *vectors)
  {
    CubicVector& cef = vectors[0];
    CubicScalar& epot = scalars[0];

    vect3d r = pos - p0;
    prec_t r2 = r * r;
    prec_t rlen = prec_t_sqrt (r2);
    if ( rlen < 1e-23 )
      {
        epot.f = q * 93 / (64 * R);
      }
    else if ( rlen >= R )
      {
        prec_t r3 = r2 * rlen;
        prec_t r5 = r3 * r2;
        prec_t r7 = r5 * r2;
        prec_t r9 = r7 * r2;

        vect3d ef = (q / r3) * r;
        prec_t f = q / r5;
        vect3d dedx = f * (vect3d(r2, 0, 0) - 3 * r.x * r);
        vect3d dedy = f * (vect3d(0, r2, 0) - 3 * r.y * r);
        vect3d dedz = f * (vect3d(0, 0, r2) - 3 * r.z * r);
        f = 3 * q / r7;
        vect3d d2edxdy = f * (5 * r.x * r.y * r - vect3d (r2 * r.y, r2 * r.x, 0));
        vect3d d2edxdz = f * (5 * r.x * r.z * r - vect3d (r2 * r.z, 0, r2 * r.x));
        vect3d d2edydz = f * (5 * r.y * r.z * r - vect3d (0, r2 * r.z, r2 * r.y));
        f = 15 * q / r9;
        vect3d d3edxdydz = f * (vect3d(r.y * r.z * r2, r.x * r.z * r2, r.x * r.y * r2) - 7 * r.x * r.y * r.z * r);

        cef.f += ef;
        cef.dfdx += dedx;
        cef.dfdy += dedy;
        cef.dfdz += dedz;
        cef.d2fdxdy += d2edxdy;
        cef.d2fdxdz += d2edxdz;
        cef.d2fdydz += d2edydz;
        cef.d3fdxdydz += d3edxdydz;

        epot.f += q / rlen;
        epot.dfdx += -ef.x;
        epot.dfdy += -ef.y;
        epot.dfdz += -ef.z;
        epot.d2fdxdy += -dedy.x;
        epot.d2fdxdz += -dedy.x;
        epot.d2fdydz += -dedz.y;
        epot.d3fdxdydz += -d2edydz.x;
      }
    else
      {
        prec_t r4 = r2 * r2;
        prec_t r6 = r4 * r2;
        prec_t r8 = r6 * r2;
        prec_t M = f35_8R5 * r2 - f21_4R7 * r4 + f15_8R9 * r6;
        vect3d ef = q * M * r;
        prec_t N = f35_4R5 - f21_R7 * r2 + f45_4R9 * r4;
        vect3d dedx = q * (r.x * N * r + vect3d(M, 0, 0));
        vect3d dedy = q * (r.y * N * r + vect3d(0, M, 0));
        vect3d dedz = q * (r.z * N * r + vect3d(0, 0, M));
        prec_t K = -f42_R7 + f45_R9 * r2;
        vect3d d2edxdy = q * (r.x * r.y * K * r + vect3d (r.y * N, r.x * N, 0));
        vect3d d2edxdz = q * (r.x * r.z * K * r + vect3d (r.z * N, 0, r.x * N));
        vect3d d2edydz = q * (r.y * r.z * K * r + vect3d (0, r.z * N, r.y * N));
        vect3d d3edxdydz = q * ((r.x * r.y * r.z * f90_R9) * r + vect3d(r.y * r.z * K, r.x * r.z * K, r.x * r.y * K));

        cef.f += ef;
        cef.dfdx += dedx;
        cef.dfdy += dedy;
        cef.dfdz += dedz;
        cef.d2fdxdy += d2edxdy;
        cef.d2fdxdz += d2edxdz;
        cef.d2fdydz += d2edydz;
        cef.d3fdxdydz += d3edxdydz;

        epot.f = -q * (f35_32R5 * r4 - f21_24R7 * r6 + f15_64R9 * r8 - f93_64R);
        epot.dfdx += -q * r.x * M;
        epot.dfdy += -q * r.y * M;
        epot.dfdz += -q * r.z * M;
      }
  }

  bool inContact (const vect3d& pos)
  {
    return false;
  }

private:
  vect3d p0;
  prec_t R;
  prec_t q;

  // M
  prec_t f35_8R5;
  prec_t f21_4R7;
  prec_t f15_8R9;
  // N
  prec_t f35_4R5;
  prec_t f21_R7;
  prec_t f45_4R9;
  // K
  prec_t f42_R7;
  prec_t f45_R9;
  // d3f
  prec_t f90_R9;
  // pot
  prec_t f35_32R5;
  prec_t f21_24R7;
  prec_t f15_64R9;
  prec_t f93_64R;
};


int main (int argc, char *argv[])
{
  //vect3d pos(0, 0, 0);
  //vect3d pos(0.001, 0.001, 0.001);
  vect3d pos(0.125, 0.125, 0.125);
  prec_t R = 0.75;

  prec_t fc_min, fc_max;
  prec_t pot_min, pot_max;
  {
    Statics statics;
    statics.addStaticElement (new myElement(pos, R, 1));
    ElectroDynamics dyn(statics);
    Scene scene(dyn);
    Screen screen(800, 800, 1.8, 1.8);
    screen.set_axis (vect3d(0, 0, 0), vect3d(1, 1, 0).normal(), vect3d(0, 0, 1));
    scene.calc (screen, Scene::CALC_EFIELD);
    fc_min = scene.get_fmin();
    fc_max = scene.get_fmax();
    fc_min /= 1.1;
    fc_max *= 1.1;
    scene.calibrate_field (fc_min, fc_max);
    scene.set_coloring (64, false);
    scene.render_map (screen);
    screen.write ("chargedist_efield_direct.bmp");
    scene.set_coloring (0, true);
    scene.render_lic (screen);
    screen.write ("chargedist_eflic_direct.bmp");
    scene.calc (screen, Scene::CALC_POTENTIAL);
    pot_min = scene.get_fmin();
    pot_max = scene.get_fmax();
    pot_min /= 1.1;
    pot_max *= 1.1;
    scene.calibrate_field (pot_min, pot_max);
    scene.set_coloring (32, false);
    scene.render_map (screen);
    screen.write ("chargedist_pot_direct.bmp");
  }

  MyBuilder builder(pos, R, 1);
  Octree tree;
  builder.build (tree);

  Statics statics;
  ElectroDynamics dyn(statics);
  statics.addStaticElement (new StaticOctree(tree, false, true, true, 1, 1));

  Scene scene(dyn);
  Screen screen(800, 800, 1.8, 1.8);
  screen.set_axis (vect3d(0, 0, 0), vect3d(1, 1, 0).normal(), vect3d(0, 0, 1));
  scene.calc (screen, Scene::CALC_EFIELD);
  scene.calibrate_field (fc_min, fc_max);
  scene.set_coloring (64, false);
  scene.render_map (screen);
  screen.write ("chargedist_efield_octree.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  screen.write ("chargedist_eflic_octree.bmp");

  scene.calc (screen, Scene::CALC_POTENTIAL);
  scene.calibrate_field (pot_min, pot_max);
  scene.set_coloring (32, false);
  scene.render_map (screen);
  screen.write ("chargedist_pot_octree.bmp");
}

#if 0
      {
        prec_t r4 = r2 * r2;
        prec_t r6 = r4 * r2;
        prec_t r8 = r6 * r2;
        prec_t M = 35 * r2 / (8 * R5) - 21 * r4 / (4 * R7) + 15 * r6 / (8 * R9);
        vect3d ef = q * M * r;
        prec_t N = 35 / (4 * R5) - 21 * r2 / R7 + 45 * r4 / (4 * R9);
        vect3d dedx = q * (r.x * N * r + vect3d(M, 0, 0));
        vect3d dedy = q * (r.y * N * r + vect3d(0, M, 0));
        vect3d dedz = q * (r.z * N * r + vect3d(0, 0, M));
        prec_t K = -42 / R7 + 45 * r2 / R9;
        vect3d d2edxdy = q * (r.x * r.y * K * r + vect3d (r.y * N, r.x * N, 0));
        vect3d d2edxdz = q * (r.x * r.z * K * r + vect3d (r.z * N, 0, r.x * N));
        vect3d d2edydz = q * (r.y * r.z * K * r + vect3d (0, r.z * N, r.y * N));
        vect3d d3edxdydz = q * ((r.x * r.y * r.z * 90 / R9) * r + vect3d(r.y * r.z * K, r.x * r.z * K, r.x * r.y * K));

        cef.f += ef;
        cef.dfdx += dedx;
        cef.dfdy += dedy;
        cef.dfdz += dedz;
        cef.d2fdxdy += d2edxdy;
        cef.d2fdxdz += d2edxdz;
        cef.d2fdydz += d2edydz;
        cef.d3fdxdydz += d3edxdydz;

        epot.f = -q * (35 * r4 / (32 * R5) - 21 * r6 / (24 * R7) + 15 * r8 / (64 * R9) - 93 / (64 * R));
        epot.dfdx += -q * r.x * M;
        epot.dfdy += -q * r.y * M;
        epot.dfdz += -q * r.z * M;
      }
#endif
